import math
from pylab import *
import time
import struct
import numpy as np
from matplotlib import pyplot

def anint(x, lx):
    lhx=0.5*lx;
    while x>lhx:
        x-=lx
    while x<-lhx:
        x+=lx
    return x

def distancia(lbox,pos1,pos2):
    dis2=0
    drx=pos2[0]-pos1[0]
    drx=anint(drx,lbox[0])
    dry=pos2[1]-pos1[1]
    dry=anint(dry,lbox[1])
    drz=pos2[2]-pos1[2]
    #drz=anint(drz,lbox[2])
    dis2=drx*drx+dry*dry+drz*drz
    return  [math.sqrt(dis2),drx,dry,drz]

def Lboundary(x, lx):
    res=x
    while x>=lx:
        x-=lx
    while x<0:
        x+=lx
    return x

def boundary(pos, lbox):
    pos[0]=Lboundary(pos[0],lbox[0])
    pos[1]=Lboundary(pos[1],lbox[1])
    #pos[2]=Lboundary(pos[2],lbox[2])
    return pos

def binaryExtract(f):
    s = f.read(8)
    t= struct.unpack("d", s)[0]
    s = f.read(4)
    nacc = struct.unpack("i", s)[0]
    s = f.read(4)
    npart = struct.unpack("i", s)[0]
    s = f.read(8*3)
    lbox= struct.unpack("ddd", s)
    posiciones=[]
    velocidades=[]
    fuerzas=[]
    diametros=[]
    #Comienza el bucle sobre el numero de particulas
    for i in range(npart):
        #Extraemos las posiciones
        s = f.read(8*3)
        pos= list(struct.unpack("ddd", s))
        #Extraemos las velocidades
        s = f.read(8*3)
        vel= struct.unpack("ddd", s)
        #Extraemos las fuerzas
        s = f.read(8*3)
        force= struct.unpack("ddd", s)
        #Extraemos el diametro
        s = f.read(8)
        diam= struct.unpack("d", s)[0]
        boundary(pos,lbox)
        posiciones.append(pos)
        velocidades.append(vel)
        fuerzas.append(force)
        diametros.append(diam)
    return [npart, nacc, lbox, posiciones,velocidades,fuerzas,diametros]

#def histogram(minValue,maxValue,values, nhist):
#    suma=0.0
#    
#    dist = np.zeros((nhist+1))
#    count = np.zeros((nhist+1))
#    step=(maxValue-minValue)/float(nhist)
#    for val in values:
#        if val>=minValue and val<=maxValue:
#            dnd=int(math.floor((val-minValue)/step))
#            dist[dnd]+=val
#            count[dnd]+=1.0
#    
#    #Normalizacion
#    for i in range(1,len(count)):
#        suma+=step*(count[i]+count[i-1])*0.5
#    x=[]
#    y=[]
#    for i in range(len(count)):
#        x.append(minValue+(float(i)+0.5)*step)
#        y.append(count[i]/suma)
#    return [x,y]

def histogram(minValue,maxValue,values, nbins):
    suma=0.0
    toRad=math.pi/180.0;
    dist = np.zeros((nbins))
    count = np.zeros((nbins))
    step=(maxValue-minValue)/float(nbins)
    print step,nbins
    for val in values:
        if val>=minValue and val<=maxValue:
            dnd=int(math.floor((val-minValue)/step))
            dist[dnd]+=val
            count[dnd]+=1.0
            suma+=1.0
    
    x=[]
    y=[]
    sss=0.0
    for i in range(len(count)):
        x.append(minValue+float(i+0.5)*step)
        fracArea=(math.cos(toRad*(x[i]-0.5*step))-(math.cos(toRad*(x[i]+0.5*step))))*0.5
        y.append(count[i]/fracArea/suma)
    return [x,y]
        
#def conectividad(rcontact, pos, lbox):
#    npart = len(pos)
#    angulos=[]
#    for i in range(npart-1):
#        for j in range(i+1,npart):
#            [dis,dx,dy,dz] = distancia(lbox,pos[i],pos[j])
#            if dis<rcontact and pos[i][2]>5.0 and pos[i][2]<lbox[2]-5.0 and pos[j][2]>5.0 and pos[j][2]<lbox[2]-5.0:
#                theta=math.acos(dz/dis)*180.0/math.pi
#                angulos.append(theta)
#                angulos.append(180.0-theta)
#    return angulos

def conectividad(pos, lbox, rcontact, zseg):
    npart = len(pos)
    angulos=[]
    #Recorremos todos los pares de particulas
    for i in range(npart-1):
        for j in range(i+1,npart):
            [dis,dx,dy,dz] = distancia(lbox,pos[i],pos[j])
            if dis < rcontact and pos[i][2]<(lbox[2]-zseg) and pos[i][2]>zseg and pos[j][2]<(lbox[2]-zseg) and pos[j][2]>zseg:
                #Las particulas estan en cntacto-> extraemos theta entre i-j
                theta=math.acos(dz/dis)*180.0/math.pi
                angulos.append(theta)
                #Por simetria incluimos angulo entre j-i
                angulos.append(180.0-theta)
    return angulos

name='/home/jcfernandez/remote/rperi76/simulaciones/eps1T0.1/u100/AMPL/W2.0/A20.0/config.bin'
name='/home/jcfernandez/remote/almacen01/jcfernandez/simulaciones/espresso/R-ShiftedLJ/phi0.2/T0.1/u20/AMPL/W2.0/A60.0/config.bin'
#name='/home/jcfernandez/remote/rperi76/simulaciones/eps1T0.1/u20/AMPL/W2.0/A0.5/config.bin'
#name='conectA100.dat'
gamma0=0.01
omega=2.0
ymax=10
rcontact=1.05
x2 = arange(0,10.0,0.01)
n=0
nnn=0
pt=None
ion()
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

zseg=1.5
f = open(name,'r')
xxx=0.5*180.0/200.0
for j in range(1000):
    [npart, nacc, lbox, posiciones,velocidades,fuerzas,diametros]=binaryExtract(f)
    
    print n
    if n==0:
        angulos = conectividad(posiciones, lbox, rcontact, zseg)
        [x,y]=histogram(0.,180.,angulos, 200)
        ax.bar(x,y)
        x1=x[:]
        y1=y[:]
        ax.plot(x1,y1,'r-')
        ax.axis([0,180,0,ymax])
    if n>600:
        
        angulos = conectividad(posiciones, lbox, rcontact, zseg)
        [x,y]=histogram(0.,180.,angulos, 200)
#        for k in range(len(y)):
#            vol=(math.pi*1.05*1.05*abs(math.cos((x[k]-xxx)*math.pi/90.)-math.cos((x[k]+xxx)*math.pi/90.)))/2.
#            y[k]=y[k]/vol
        ax.clear()
        ax.bar(x,y)
        ax.plot(x1,y1,'r-')
        ax.plot([0,180],[1,1],'g-')
        ax.axis([0,180,0,ymax])
        draw()
        nnn+=1
            
    else:
        pass
#        ax.clear()
#        ax.bar(x,y)
#        ax.plot(x1,y1,'r-')
#        ax.axis([0,180,0,ymax])
#        draw()
    n+=1
f.close()
